function [DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeMan(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,natio,singl)
global dimS dimProv dimNum
% auxG=[GAMMA';1-GAMMA'];
% auxG=kron(auxG,ones(dimS,1));
% auxG=reshape(auxG,2*dimS*dimProv*dimNum,1);
auxG=[repmat(GAMMA'.*kron(Weight,ones(1,dimNum)),[dimS 1]);repmat((1-GAMMA').*kron(Weight,ones(1,dimNum)),[dimS 1])];
auxG=auxG(:);

auxSurp0=[STORE;STORE_];
auxSurp0=permute(auxSurp0,[1 3 2]);auxSurp0=reshape(auxSurp0,2*dimS*dimProv*dimNum,14);
if singl==-1;
    if natio>0
    ind=(auxSurp0(:,1)==natio);
    elseif natio<0
      ind=(auxSurp0(:,1)~=abs(natio));  
    elseif natio==0
      ind=logical(ones(2*dimS*dimProv*dimNum,1));  
    end
else
    if natio>0
    ind=(auxSurp0(:,1)==natio)&(auxSurp0(:,7)==singl);
    elseif natio<0
    ind=(auxSurp0(:,1)~=1)&(auxSurp0(:,1)~=abs(natio))&(auxSurp0(:,7)==singl);
    elseif natio==0
    ind=logical(ones(2*dimS*dimProv*dimNum,1));  
    end
end
auxSurp1=[STORE1;STORE1_];
auxSurp1=permute(auxSurp1,[1 3 2]);auxSurp1=reshape(auxSurp1,2*dimS*dimProv*dimNum,14);
e=(1:dimProv*dimNum*dimS*2)';
Identity0=e(ind);
ss0=zeros(length(Identity0),1);
ss1=zeros(length(Identity0),1);
ss0_L=zeros(length(Identity0),1);
ss0_H=zeros(length(Identity0),1);
ss1_L=zeros(length(Identity0),1);
ss1_H=zeros(length(Identity0),1);
Drho=zeros(length(Identity0),1);
DS1=zeros(length(Identity0),1);
DS2=zeros(length(Identity0),1);
DS3=zeros(length(Identity0),1);
DS4=zeros(length(Identity0),1);

ww=zeros(length(Identity0),1);
parfor i=1:length(Identity0)
    id0=Identity0(i);
    id1=id0;
    ss0_L(i)=auxSurp0(id0,9);  % Lowest utility in marriage is outside option
    ss0_H(i)=(1-auxSurp0(id0,7)).*(auxSurp0(id0,11)-auxSurp0(id0,10))+auxSurp0(id0,7).*auxSurp0(id0,9);
    ss0(i)=(1-auxSurp0(id0,7)).*auxSurp0(id0,11).*auxSurp0(id0,9)./(auxSurp0(id0,9)+auxSurp0(id0,10))+auxSurp0(id0,7).*auxSurp0(id0,9);
    ss1_L(i)=auxSurp1(id1,9);
    ss1_H(i)=(1-auxSurp1(id1,7)).*(auxSurp1(id1,11)-auxSurp1(id1,10))+auxSurp1(id1,7).*auxSurp1(id1,9);
    ss1(i)=(1-auxSurp1(id1,7)).*auxSurp1(id1,11).*auxSurp1(id1,9)./(auxSurp1(id1,9)+auxSurp1(id1,10))+auxSurp1(id1,7).*auxSurp1(id1,9);
    Drho(i)=auxSurp1(id1,9)./(auxSurp1(id1,9)+auxSurp1(id1,10))-auxSurp0(id0,9)./(auxSurp0(id0,9)+auxSurp0(id0,10)); % change in the weight
    DS1(i,1)=(1-auxSurp0(id0,7)).*(1-auxSurp1(id1,7)).*auxSurp1(id1,11)*Drho(i); % Stay married, change in weight
    DS2(i,1)=(1-auxSurp0(id0,7)).*(1-auxSurp1(id1,7))*auxSurp0(id0,9)./(auxSurp0(id0,9)+auxSurp0(id0,10))*(auxSurp1(id1,11)-auxSurp0(id0,11)); % D surplus
    DS3(i,1)=auxSurp0(id0,7).*(1-auxSurp1(id1,7))*(ss1(i)-auxSurp0(id0,9)); % Single to married
    DS4(i,1)=(1-auxSurp0(id0,7)).*auxSurp1(id1,7)*(auxSurp1(id1,9)-ss0(i)); % Married to single
    ww(i)=auxG(i);
end
ww=ww/sum(ww);
S0=[ss0_L ss0_H ss0_L ss0_H];
S1=[ss1_L ss1_L ss1_H ss1_H];
DSurplus=S1-S0;
[DSurplusMin,i2]=min(DSurplus,[],2);
i3=sub2ind(size(DSurplus),(1:length(Identity0))',i2);
[DSurplusMax,i2]=max(DSurplus,[],2);
i4=sub2ind(size(DSurplus),(1:length(Identity0))',i2);

DSurplus_L=sum(DSurplusMin./S0(i3).*ww);
DSurplus_H=sum(DSurplusMax./S0(i4).*ww);
DSurplus=sum((ss1-ss0)./ss0.*ww);
Oaxaca=sum(ww.*[DS1 DS2 DS3 DS4]./ss0);
%Oaxaca=Oaxaca/sum(Oaxaca)*100;